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We provide a simple but robust argument that primordial black hole (PBH) production generically 
does not exceed astrophysical bounds during the resonant preheating phase after inflation. This 
conclusion is supported by fully nonlinear lattice simulations of various models in two and three 
dimensions which include rescattering but neglect metric perturbations. We examine the degree to 
which preheating amplifies density perturbations at the Hubble scale and show that at the end of 
the parametric resonance, power spectra are universal, with no memory of the power spectrum at 
the end of inflation. In addition we show how the probability distribution of density perturbations 
changes from exponential on very small scales to Gaussian when smoothed over the Hubble scale 
- the crucial length for studies of primordial black hole formation - hence justifying the standard 
assumption of Gaussianity. 
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I. INTRODUCTION 

Primordial black holes (PBHs) span a wide range of 
mass scales and are typically much smaller than the so- 
lar mass (~ 10 33 g) and may be formed in the early uni- 
verse . PBHs may form from the gravitational collapse 
of large density fluctuations at horizon (i.e. Hubble scale 
k = aH) crossing in the radiation dominated universe. 
A PBH formed at the Planck time ~ 10 _43 sec will have 
a mass ~ 10 19 GeV, while masses around ~ 10 15 g are 
formed at ~ 10 _23 sec ( see for example Q ). 

The evaporation time for a PBH mass ~ 10 15 g is nearly 
the present age of the universe, so PBH with smaller 
masses than this would have evaporated in the past, un- 
loading a potentially vast amount of entropy. The success 
of the standard cosmology and the observation of the cos- 
mic rays severely constrains the abundance of PBHs for 
various masses and provides useful constraints on infla- 
tionary and early universe physics. 

For example, Big Bang Nucleosynthesis limits the PBH 
abundance in the mass range 10 6 ~ 10 13 g, by limiting 
the entropy from Hawking radiation or requiring that it 
does not modify the cosmological composition of the light 
elements Q- PBH of mass 10 15 g evaporating now will 
emit particles such as 7-rays which are constrained by 
observations of the extragalactic 7-ray background which 
imply VlpBH < 1CT 8 0. For > 10 15 g, the limit on the 
PBH abundance is obtained from requiring that £Ipbh 
does not exceed unity. 

These observational constraints on PBHs provide a 
powerful probe of the primordial fluctuations. The upper 
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bound on the abundance of PBHs directly leads to that 
on the density fluctuation at horizon crossing when PBH 
are formed. Therefore the scales of fluctuations relevant 
to PBH formation are much smaller than those associ- 
ated with the Cosmic Microwave Background and the 
large-scale structure. This makes studying PBHs impor- 
tant. In the past, for example, constraints on the density 
perturbation spectrum were obtained by studying PBH 
formation While the requirement that PBHs are not 
over-produced yields useful information about the early 
universe, PBH can be an interesting dark matter candi- 
date in smaller abundances @- 

In this paper we will show that typically PBH are not 
over-produced during the violent non-equilibrium phase 
of preheating that follows the end of many inflationary 
models. This follows from three key observations: (1) 
the peak of the density perturbation spectrum typically 
lies at scales smaller than the Hubble scale. (2) The peak 
corresponds to density contrasts of order unity. (3) The 
slope of the spectrum around the horizon size is three 
(in 3d). Putting these together we typically find that at 
the horizon scale relevant for PBH formation, the den- 
sity contrast is around an order of magnitude too small 
to over-produce PBH. Nevertheless, the density pertur- 
bation on the horizon scale is significantly enhanced by 
preheating (by several orders of magnitude) compared 
with the no-resonance case and hence preheating is im- 
portant in understanding the potential astrophysical and 
cosmological implications of PBH. 

In studying the production of PBHs, one usually as- 
sumes that the probability distribution of density fluctu- 
ations at horizon crossing is Gaussian. This assumption 
is critical because the density perturbations which col- 
lapse to black holes are very rare: several a fluctuations 
(otherwise PBH will be over-produced in all cases) and 
therefore production of PBHs is sensitive to the tail of 
the distribution. Indeed Bullock and Primack found 
that in some inflation models large perturbations are sup- 
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pressed relative to a Gaussian distribution, resulting in 
a significant change in a number of PBH. We study the 
validity of the Gaussian assumption in a later section. 

Among various possible scenarios that might over- 
produce PBHs, we focus on preheating after inflation. 
Preheating is a process in which energy transfer occurs 
rapidly from inflaton field to another field due to the non- 
perturbative effects during the oscillating phase of the 
inflaton This process significantly differs from the 

usual reheating scenario, where inflaton decays perturba- 
tively to another particles, in a sense that in preheating 
most energies of inflaton field converts to created parti- 
cles only during the several oscillations of inflaton and 
even the massive particle which is much heavier than the 
inflaton can be created. It has been understood that 
parametric resonance occurs generically at the first stage 
of reheating ^l|. The parametric resonance does not 
last long because the rapid increase of created particles 
eventually affect the motion of background field and the 
created particles scatter off each other, removing the par- 
ticles from the resonance band. By these effects, the res- 
onance becomes inefficient and the decaying process of 
inflaton is described by usual single-body decay theory 
which finally leads to the thermal equilibrium state [l2fl . 

There are many works on preheating. The first stage 
of preheating, where the backreaction is negligible and 
the linear approximation is valid, was studied in detail 
in |l3t Il4l Il6| . Parametric resonance including metric 
perturbations in order to see the behavior on super hori- 
zon scales were studied in [IE E3, E3, 113, ISE H^- 
As in the case without metric perturbations, there is 
a crucial difference between the single and multiple 
field cases. The analysis of fully non-linear preheat- 
ing including gravity is very difficult. Before now, dif- 
ferent approximations such as mean-field approxima- 
tion |M El IH IH HHHm, lattice simulation HJ 
without metric perturbation and one dimensional fully 
non-linear calculations |3lL I32I l33| have been done for 
studying the various effects caused by preheating on the 
present universe. 

Green and Malik |34[ argued that PBHs will be over- 
produced due to the amplification of fluctuations dur- 
ing preheating for many parameter regions in a two-field 
massive inflation model based on the results of j3f| , which 
takes into account of the second order fluctuation of \ 
field. Put simply, their results suggested that the back- 
reaction timescale was smaller than the timescale for the 
over-production of PBH. 

Bassett and Tsujikawa [2(j studied PBH production 
in the two-field massless inflation model including the ef- 
fect of backreaction via the Hartree-Fock approximation 
and found that PBH overproduction might occur, if the 
probability distribution of fluctuation at horizon crossing 
was assumed to be chi-squared (which lowers the thresh- 
old mass variance, a). However, they found that PBH 
were not over-produced if the distribution was assumed 
to be Gaussian and the density field was smoothed on the 
scale of the horizon. Nevertheless, that analysis was lim- 



ited since it neglected the mode-mode coupling effects of 
rescattering and hence there was an open question both 
as to the underlying probability distribution of the den- 
sity fluctuations and the contribution of rescattering to 
horizon-scale fluctuations. 

We address both of these issues in this paper. We 
study PBH production due to preheating via two and 
three-dimensional lattice simulations which automati- 
cally include the effects of backreaction and rescatterings. 
We modified the C++ code LATTICEEASY written by 
Felder and Trachev j3||. In lattice simulations, the evo- 
lution equations for the scalar fields (and also the scale 
factor) are solved in real (as opposed to Fourier) space 
(N = 2048 2 for 2 dimensional simulation and N = 128 3 
for 3 dimensional simulation) . Metric perturbations were 
not included, so we cannot apply this method to the dy- 
namics on super horizon scales. We followed the evolu- 
tion of both the scalar fields and the total density pertur- 
bation. We found that PBH are not overproduced and, 
interestingly, that the power spectrum at the end of pre- 
heating has a universal feature, that is, it is determined 
by the preheating dynamics and does not depend on the 
initial conditions. We also studied the probability dis- 
tribution of the density perturbation at horizon crossing 
and found that it remained Gaussian which seems to be 
valid even in the tail of distribution. 

A brief summary of our paper is as follows. Section 
II gives a brief review of calculating the PBH abundance 
formed from the large density perturbation. Section III 
describes the model of preheating in this simulation. In 
section IV, we consider the initial conditions of scalar 
fields. Section V shows the numerical results of lattice 
simulation. We will interpret the non-linear behavior of 
fluctuations and study the production of PBHs. Section 
VI discusses the probability distribution of fluctuations 
at horizon crossing we assume Gaussian in usual and Sec- 
tion VII is a conclusion. 



II. PBH FORMATION BY PARAMETRIC 
RESONANCE 

A. Abundance of PBHs 

In this subsection, we briefly review the standard 
method to estimate the abundance of PBHs [3]} • In the 
radiation dominated universe PBHs will be produced if 
Sh , the amplitude of density perturbations 5 smoothed 
over the horizon size in the comoving gauge, exceeds a 
certain threshold 5 C |3Sl l39j . From linear analysis the 
critical value of <5 C is roughly estimated to be 1/3, but it 
is not independent of the initial density profile. Numeri- 
cal study |4(| suggest 5 C ~ 0.7 for various initial density 
profiles (see also |4j|). 

Under the assumption that the probability distribution 
of density fluctuation at horizon crossing is Gaussian, 
the mass fraction of PBHs at the formation time, /3, is 
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estimated as 

/oo / ^2 \ 

,,**'^"(^)' "> 

where a 2 is the variance of 8h- Since the mass fraction 
of PBHs increases in proportion to the scale factor in the 
radiation dominated universe, (3 must be very small in 
order to satisfy the astrophysical constraints. 

Roughly speaking, /3 is observationally constrained to 
be smaller than about 10 -20 for most of the range of PBH 
mass except for a small window at M » 10 15 g, which 
corresponds to PBHs evaporating now. If we adopt 10~ 20 
as the upper bound on (3, the upper bound on a becomes 
~ 0.03 and ~ 0.08 for S c — 0.3 and 0.7, respectively, 
assuming a Gaussian distribution for S. Though there is 
uncertainty in 8 C it does not affect our conclusions in the 
range of values 0.3 — 0.7. Therefore we adopt the smaller 
value 0.03 as the upper bound on a to be conservative. 

To estimate the variance of density perturbations at 
the horizon scale, we simply use the power spectrum 
Vs (aH) , where a and H are, respectively, the scale factor 
and the Hubble parameter, and Vs(k) is defined by 

(6 k 6^) = ^-V s (k)S(k-k'). (2) 



B. Models of parametric resonance 

In this paper, we consider two simple models of pre- 
heating. 



1. Conformal Models 

Conformal models are models composed of two scalar 
field with the potential given by [lj| , 

v(<t>,x) = 2<f + (3) 

Here is the inflaton field. We start our simulation 
at the time when drops down to 0.34m p / |l4j . where 
to 2 ; = 1/G. In the oscillating phase of the inflaton the 
universe is effectively radiation dominated for the poten- 
tial quartic in fields when averaged in time. 
As standard, we introduce rescaled fields by 

= a^/0o, x = a x/0o, (4) 

where the scale factor a is normalized to unity at the end 
of inflation, (i.e., at the beginning of preheating) , 

We also introduce the rescaled conformal time fj, re- 
lated to proper time t by 



Then, the equations of motion for this model become 

0"-A- 1 A0 + 3 + ^0x 2 --0 = O, (6) 
A a 

x"-A- 1 Ax+^0 2 x--X-O, (7) 
A a 

where ' denotes the differentiation with respect to fj. In 
the radiation dominated universe ( a oc fj) the last terms 
in Eqs. JHJl and Q proportional to a" vanish, and thus 
these equation reduce to the Minkowski ones. This is 
only exactly true if were conformally coupled to the 
curvature but since this is a weak effect we neglect it. 

The unperturbed background solution for is given 
by Jacobi's elliptic cosine function, 01(77) Then, lin- 
earized equations obey the so-called Lame equation |43| 
with resonance parameters 3 and g 2 /X for and Xi re ~ 
spectively. Hence, the growth rate of the longest wave- 
length mode for \ is solely determined by g 2 /X, and 
there is a strong resonance at the longest wavelengths 
for g 2 /X — 2,8,18, Roughly speaking, the largest 
wave number in the efficient resonant band is 

/ g 2 \ i/4 

k max = ( y J VA0O, (8) 

An outstanding feature of conformal models is that the 
modes which are amplified by parametric resonance do 
not change by cosmic expansion. The Lame equation 
for also has instability bands, but the growth rate is 
small compared with the typical one for x and limited to 
roughly the Hubble scale 0] . 

In our lattice simulations the evolution of the scale fac- 
tor a is determined self-consistently by solving Friedmann 
equation with the spatially averaged energy density. In 
the simulation, A is fixed to 9 x 10~ 14 appropriate to the 
COBE normalization. We studied both g 2 /A = 2 and 
50. In both cases there is parametric resonance of x f° r 
k = mode in the linear regime. In the former, the x 
background is not strongly suppressed while in the lat- 
ter case the x ne ld is heavy and is strongly suppressed 
during inflation [21121 12^ . 

2. Massive Inflaton Models 

We also considered massive inflaton models with the 
potential 

2 2 
Tfl Q 

V(<t>,x) = -2 2 + y0V- (9) 

When the inflaton field oscillates around the potential 
minimum, the equation of state of the inflaton is dust on 
average. This means that amplitude of the background 
inflaton field decreases as oc a~ 3//2 . Therefore it is con- 
venient to introduce rescaled fields as 



adfj = "vX<podt. 



(5) 



= a 3 / 2 0/0o, X = a 3/2 XM), (10) 
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where </>o = 0.193m p /. Then, until the back reaction due 
to parametric resonance becomes efficient, the amplitude 
of the background <j> stays almost constant. In terms of 
rescaled fields the equations of motion are 

4> - a- 2 A4> - ^H<j> - ^H 2 4> + m 2 + g 2 a- 3 4>x 2 = 0, 

X - a- 2 Ax - ^Hx ~ \h 2 X + ,g 2 a^ 2 X = 0. 

(11) 

where ' denotes the differentiation with respect to the 
proper time. From these equations, in contrast to the 
conformal case, we see that the expansion of the universe 
affects the motion of scalar fields: the wavelength of each 
comoving mode is redshifted and the effective coupling 
between <fi and x 1S decreased. 

The linearized equations of Eqs. 1111) was extensively 
studied in |l3j| . The equation for x approximately reduces 
to the so-called Mathieu equation and the evolution of x 
field shows broad resonance for 



q 



9^0 
4m 2 



> 1. 



(12) 



Huge amplification of x occurs at each time when ampli- 
tude of the inflaton field becomes zero due to violation 
of the adiabatic condition lu/lo 2 < 1 for the x-field fre- 
quency, u>. 

When the expansion of the universe is taken into ac- 
count the parametric resonance shows stochastic behav- 
ior because the phase of the x ne ld is randomized. De- 
spite the stochastic nature, the amplitude of x field grows 
exponentially on average. The maximum wave number 
in the resonance band is given bv|l3fl. 



1/4 
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In our simulations we adopt m = 10 



(13) 

,i indicated by 
the COBE normalization, and g = 1CP 3 as a represen- 
tative which realizes strong resonance. With this choice 
of parameters, we have q = 10 at the beginning of pre- 
heating. 



is the effective mass evaluated at the end of inflation. 
Recall that we have set a = 1 at that time. Since we 
consider the cases with m x /k max — (g 2 /X) 1 ^ 4 greater 
than unity ( for g 2 /X = 50), the power spectrum of x 
field is given by 



V x {k) 



(16) 



for wavelengths relevant for parametric resonance. 

We also use Eq. (|rl|) as the initial condition for fluctu- 
ations of </), replacing m x in tOk with 



m<f, 



3A0o, 



(17) 



in the conformal models and m in the massive inflation 
models. 

In the conformal model with g 2 /\ = 2, x stays almost 
always massless during inflation. In this case the initial 
spectrum on the superhorizon size becomes scale invari- 
ant |23ll25j |. However the precise initial power spectrum ( 
i.e. after inflation) of scalar fields on subhorizon scales is 
rather involved. So here, for simplicity, we approximated 



V 



(18) 



which is blue on small scales and flat on large scales, 
capturing the key features of the spectrum. We will see 
that the precise form has no effect on the final results. 

In our simulations, we compute classical dynamics tak- 
ing the variance of these initial quantum fluctuations as 
if it were statistical variance, as standard [36| . Dur- 
ing the early stage the evolution is in linear regime and 
amplification of the x occupation number in the quan- 
tum picture is correctly described by the amplification 
of the perturbation amplitude in classical dynamics |l7j . 
When the nonlinearity becomes important, the occupa- 
tion number of modes relevant for resonance is far be- 
yond unity. Hence, a classical treatment is justified at 
late epoch, too. Moreover as we will see later, our final 
conclusion is quite insensitive to initial conditions. 



C. Initial Conditions 

Since the modes that our simulation covers are mostly 
on subhorizon scales, the initial conditions for the field 
fluctuations after inflation can be determined by the for- 
mula for the adiabatic vacuum. For example, for the 
X-field we have 

(XkX-v) = ^-S(k - k'), (14) 



where u>k = y m 2 , + k 2 and 

m x := # , (15) 



III. NUMERICAL RESULTS 
A. Conformal Models 

1. <? 2 /A = 50 case 

Fig. n shows the time evolution of the power spectra 
of the </> and x fields and density perturbations until the 
backreaction shuts off the parametric resonance. Let us 
first focus on the power spectra of and x fields. Both 
the power spectra are proportional to k 3 at the beginning. 
This is because the mass of each field is larger than the 
highest momentum resolved by simulation. In the early 
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FIG. 1: Evolution of the power spectrum of scalar fields and 
the total density perturbation for g 2 /A = 50. The horizontal 
axis is comoving wave number normalized by the (effective) 
inflaton mass at the end of inflation m0. Time flows ver- 
tically. Lines labeled as a, b, c, d, e and f are for times 
fj = 0, 50, 60, 70, 80, 100 respectively. Thus during our range 
of simulation, the inflaton field oscillates about 100/7.4 ~ 14 
times. An arrow pointing upward represents rapid increase of 
the amplitude of fluctuations of the inflaton field due to the 
initial parametric resonance, and the arrow pointing right- 
ward shows rescattering: high energy particles are generated 
by mode-mode coupling and scattering of low energy parti- 
cles. 



stage of evolution linear perturbation is a good approxi- 
mation. The maximal characteristic (Floquet) exponent 
Umax for <p is ~ 0.036, while that for x is ~ 0.2. Hence, 
in this early stage perturbations of ^-field grow exponen- 
tially, but those of </>- field almost stay constant. After 
a few oscillations of the inflaton field, the perturbations 
of (j> suddenly start to grow. This can be understood 
as follows. From Eq. 10 , the equation of motion for 
(4>) is 

(19) 

As x grows exponentially by the parametric resonance, 
the last term in Eq. ljT$|) . ^(4>}x , exceeds the term 

~ 2 

3(</>) (p. At this stage, the linear approximation for (p 
breaks and (p begins to grow proportional to x 2 , whose 
characteristic exponent is ~ 0.4 [lj|. This happens when 
X exceeds a critical value, 



Xc ~ rr 



-1/2 



(20) 



At the largest wave number in the resonance band 
of ^-field, the initial amplitude of (p is given by » 

{g 2 /\f/ & s/\f^. Hence Xc is esti- 



£■3 



"max/^VX 

mated as 



-1/16 



A 1 / 4 ~ 10" 



(21) 



This rough estimate is consistent with the results of our 
lattice simulations. 

Exponential amplification due to parametric resonance 
still continues until the effect of backreaction becomes 
significant. Parametric resonance ends when the second 
term in Eq. (JSJ becomes equal to the first term, that is 
when x is amplified to 



Xb 



-1/4 



(22) 



The initial amplitude of x at the shortest resonant mode 
given in © is » y/k max /m x w (g 2 /X) 1 ^ 8 y/\. Approxi- 
mating the evolution of % as x eM '' > the time at which 
parametric resonance ends, 77/, is estimated by 

-3/8 



A 



(23) 



Solving Eq. (|23|l . we have fjf ~ 70. This estimate is con- 
sistent with the result ijf ~ 80 read from Fig.^ Since fjf 
depends on the initial amplitude of x logarithmically, fjf 
does not depend so much on the parameter g 2 /X as long 
as models associated with strong resonance are concerned 

m 

By the time when the backreaction becomes impor- 
tant, the amplitude on smaller scales also increases due 
to the effect of rescattering. 
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In three dimensions, at the time when the simulation 
ends, the shortest wavelength modes have the largest am- 
plitude in the simulation. However, this seems to be an 
artifact due to lack of resolution. In the corresponding 
two dimensional simulations shown in FigEl we can also 
see the rescattering effect. In this case perturbations do 
not pile up near the shortest wavelength but peak at a fi- 
nite value of k. This is indeed consistent with the picture 
that x particles with very large kinetic energy k — > oo, 
are not produced. 

Let us now focus on the power spectrum of density 
perturbations. The energy density p is given by 

« 4 a-V = - nif + - Hxf + ^(v*) 2 

where TL := log a. From this equation, the density 
perturbations to first order are 

ctX-Hp a {{$)' - H$)){tft - fbp) + (0) V (25) 

As we have already mentioned, amplitude of tp is almost 
constant in the linear perturbation regime. Therefore 
in this regime amplitude of density perturbations does 
not grow. This behavior of Sp can be observed in Fig.^ 
When the amplitude of \ reaches x Cl the terms second or- 
der in field perturbations start to contribute to Sp. Then, 
density perturbations begin to grow rapidly. After the 
exponential increase of Sp, the amplification of density 
perturbations stops when the growth of field perturba- 
tions terminates due to back reaction to the oscillation 

of (4>). 

From Fig. ^ we find that the slopes of power spectrum 
of resultant density perturbations on large scales are all 
equal to three. 

As we shall see below, this result does not change even 
if we artificially amplify the initial fluctuations of fields 
on large scales, which leads to the rule of thumb that af- 
ter preheating correlation of perturbations on large scales 
disappears rather independently of initial power spec- 
trum. 

We can give a simple interpretation to this result. 
What we assume is that the resultant density pertur- 
bations have a typical scale r, and correlations on larger 
length scales are strongly suppressed. In such a situation, 
the integral 

^jVsi^oc J d D x(5(0)5(x))e^, (26) 

is dominated by a small region with |x| < r, where D 
is the number of spatial dimensions, which is 3 in our 
simulation. For long wavelength modes with k <C r^ 1 , 
e lkx can be approximated by unity. Thus we have 

Vs(k) cx k D , (27) 

We also performed 2 dimensional lattice simulations, and 
found the power spectrum proportional to k 2 on large 
scales as predicted; see (Fig. [2J • 
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FIG. 2: The power spectrum of density perturbation after 
preheating for two-dimensional space. We see that the slope 
of the power spectrum for small k is two, as expected. 

One may think that we have obtained this result be- 
cause of the initial blue spectrum of x-field. Since the 
characteristic exponent of the parametric amplification 
is almost the same for all modes with kjm^ < (g 2 /X) 1 ^ 4 , 
the parametric resonance will end due to the backreaction 
from the mode of the shortest resonance scale, at which 
the initial amplitude is the largest among the modes in 
resonance. At that time perturbations of x-filed still re- 
main small on large scales (where by large scale we here 
mean around the horizon size and larger). 

Here we show that the initial blue spectrum is not a 
necessary condition for suppression on large scales. For 
this purpose, we performed the same simulation but with 
the scale invariant initial spectrum (7 , x =constant, where 
the power spectrum is initially amplified on large scales) . 
Fig. U| shows Vs after preheating in this case. 

From this figure, we see that at the end of preheating 
perturbations on horizon scales are suppressed with slope 
3, which is the same as in the case with the initial blue 
spectrum. This result indicates that in general paramet- 
ric resonance causes loss of correlation between density 
fluctuations beyond a typical length 

scale and energy is efficiently cascaded to shorter wave- 
lengths by rescattering. 

In order to estimate the production rate of PBHs, we 
have to compute Sp/ p smoothed over the horizon size. 
The horizon size when parametric resonance ends corre- 
sponds to k = aH » 5 x lO~ 3 m0, which is not covered 
in our three-dimensional simulations. 

However, since there is no typical length scale before 
the horizon scale, it will be natural to expect that one can 
extrapolate the power spectrum to horizon size assuming 
the slope of the power spectrum of density perturbations 
is D. The result of 2 dimensional simulations (Fig. |2J) 
also support this extrapolation. With the aid of this 
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FIG. 3: Power spectrum of the total density perturbation 
after preheating for three cases in three dimensions. Lines 
labeled as a, b and c correspond to g 2 /X = 50, g 2 /X = 50 with 
scale-invariant power spectrum and g 2 /X = 2 respectively. 
We see that the slope of the power spectrum for small k is 
universal with a value of three. The dashed line has k 3 power 
and crosses the threshold at the horizon scale. All three lines 
are lie significantly under the dashed line, showing that PBHs 
are not overproduced in these cases. The peak of the spectrum 
and subsquent decay are not resolved by our three dimensional 
simulations but are resolved in two dimensions. 

extrapolation, the amplitude of density contrast at the 
horizon size Sh can be estimated from Fig. ^ as 

fc~3x 1(T 4 , (28) 

which is an order of magnitude smaller than the threshold 
for PBH overproduction (0.03). Therefore we conclude it 
is unlikely that PBHs arc overproduced by the parametric 
resonance in the case with g 2 /X = 50. 

2. g 2 /X = 2 case 

We also performed lattice simulations for g 2 /X = 2. In 
this case, the power spectrum of x at the end of inflation 
is fiat for k < Hq, where Ho is the Hubble parameter at 
the end of inflation. Therefore, the modes whose wave 
length is much larger than the horizon size is not sup- 
pressed in this case, which is different from the case 1 
g 2 /X 1 j2^. Taking into account the fact that there 
is a strong resonance at small k for g 2 /X = 2 , there is a 
possibility of PBHs overproduction in this case j2(| |^| . 

Vs for g 2 /X = 2 is shown in Fig. [3J As before, the 
power spectrum of i5 after parametric resonance is sup- 
pressed and its slope is three around the horizon scale, 
implying that large scale perturbations are uncorrelated. 
Therefore in the same way as discussed in the case with 
g 2 /X = 50, production of PBHs caused by parametric 



resonance will not be efficient enough to exceed the as- 
trophysical bounds. 

B. Massive Inflaton Models 

Fig. 0] shows the evolution of the power spectrum dur- 
ing preheating for a massive inflaton model. From this 
figure we find that the power spectrum after preheating is 
oc fc 3 on horizon scales as in the case of conformal models. 
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FIG. 4: Evolution of the power spectrum of S for the massive 
inflaton model. Lines labeld as a, b, c, d, e and f correspond 
to rat = 0, 50, 60, 70, 100, 120 respectively. The dashed line 
represents the threshold for PBH over-production which lies 
well above any of the curves. 

Here we cannot directly use the criterion for the PBH 
production discussed in section II, because the universe 
is not radiation dominated but dust on average. Due 
to the difference of the equation of state, the condition 
that density perturbations collapse to form a black hole 
differs from the one in the radiation dominated universe. 
In [33, the critical density S c for the equation of state 

P = ep, (29) 

depends on e. Here we assume that the instability of 
the density perturbation of scalar fields in massive in- 
flaton models is similar to the fluid case with the same 
effective equation of state. 

The energy density and pressure in massive inflaton 
models are 

p = T + U, P = T-U, (30) 

with 

T = -<j> 2 + -x 2 , 

2 2 

u = J(-v0) 2 + |(-Vx) 2 + ^ 2 + ^V.(3i) 

la la I I 
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Using the equations of motion for the scalar fields, we 
can show the relation, 

(T) = (U) + ^(^ X 2 ), (32) 

where (• • •) denotes a long time average with the weight 
a 3 dt. Hence we can estimate e by 

e«yfoV>/<P>. (33) 

Fig. shows the time evolution of this quantity during 
preheating. 




40 60 80 100 

mt 



FIG. 5: Time evolution of the ratio of the interaction energy 
4j-</> 2 X 2 to the total energy density p. Time is normalized in 
1/m. We see that e after preheating is about 3 x 10~ 2 , i.e. 
the total system behaves approximated as dust. 

At the end of preheating e becomes as large as ~ 3 x 
10~ 2 . Therefore the ratio of pressure to energy density 
after preheating in massive inflaton model is ten times 
smaller than that of the radiation dominated universe. 
Hence, the upper limit on a in the case of massive inflaton 
model will be reduced to about 3 x 10~ 3 . On the other 
hand, from Fig. 01 the value of the power spectrum of 
the density perturbation at horizon size (k w 1 x 10 m) 
can be estimated as 

S H ~ 5 x 1(T 4 , (34) 

where we have used the simulations to conclude that the 
power spectrum is proportional to fc 3 for small k. This is 
smaller than the threshold ~ 3 x 10 -3 again by about one 
order of magnitude. Therefore we conclude that PBHs 
will not be overproduced in massive inflaton model, ei- 
ther, despite the fact that a is significantly enhanced by 
resonance on horizon scales. 







^ -2 
O " 3 
-4 
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FIG. 6: Time evolution of y/Vs at fc = k max — q x ^m. From 
this, we see that backreaction terminates the parametric am- 
plification at around mt — 70. 



abundance of PBHs by assuming that the probability dis- 
tribution of density perturbations at the horizon size is 
Gaussian. If the tale of distribution, which is relevant for 
PBH formation, had a non-Gaussian tail, the resultant 
astrophysical constraints would be significantly altered 
and hence it is a crucial assumption to test. Further, 
since rescattering (oc 5\ 2 ) is crucial, one might expect 
chi-squared corrections to be important. 




FIG. 7: The probability distribution of density fluctuations 
smoothed over the scale (the Compton wavelength of the 
inflaton) in the conformal model. The dotted line is the spec- 
trum at initial time and black line is the spectrum at the 
end of preheating. The dashed line is an appropriately scaled 
Gaussian distribution. At the end of preheating, the distribu- 
tion of large fluctuations is significantly more amplified than 
the Gaussian prediction. 



IV. GAUSSIANITY 

In this section, we discuss the probablity distribution 
of the amplitude of density perturbations at the end of 
preheating. In the preceding sections we estimated the 



We first show the distribution of density fluctuation 
smoothed over the shortest resonant scale in a conformal 
model. The result is shown in Fig. [7| We can see that 
the distribution does not trace a Gaussian distribution at 
the end of preheating, while it does at the initial stage. 
In particular, the probability of large amplitude of per- 
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turbations is enhanced through preheating. Interestingly, 
the late time distribution looks like exponential. 

This results can be understood as follows. At initial 
stage where the linear approximation is valid, density 
perturbation is just a superposition of Gaussian distri- 
butions. Hence, the probability distribution is Gaussian. 
As perturbations grow, the terms quadratic in field per- 
turbations start to contribute to p. The probablity distri- 
bution in such situation will be mimiced by a product of 
two Gassian random variables x and y. The probability 
distribution of z = xy is given by 

P(z) = { 

27R71C72 J X 

1 s_ 

~ — e "I's , (35) 

2v27rcri(72Z 

where a x and <Ty are variances x and y, respectively. Here 
in the last step we used the steepest decent method as- 
suming z is much larger than o\Oi- In this manner one 
can reproduce a pure exponential distribution for large 
values of z. 




FIG. 8: The probability distribution of density fluctuation 
smoothed over the scale 50 x m$ in the conformal model. 
The dotted line is the distribution at the initial time and the 
black one is at the end of preheating. Contrary to Fig. [7] the 
distribution is still Gaussian even at the end of preheating. 
Hence, on the horizon scale relevant to PBH production the 
Gaussian assumption is a good approximation. 

Next we show the distribution of density perturbations 
averaged over a large scale. The result is shown in Fig. [8] 
In this case the distribution is almost Gaussian even at 
the end of preheating. This result can be interpreted 
as follows. As we discussed in section V, the density 
perturbations lose correlation on scales much larger than 
the shortest wavelength in the resonance band. Hence, 
the average over a large length scale L behaves as a sum 
of a large number of independent random variables of 
O ((Lkmax/a) 3 )- Therefore its distribution is guaranteed 
to be close to Gaussian by the central limit theorem, 
which is consistent with numerical the results. Significant 



deviations from Gaussianity are not expected unless the 
amplitude is as about (Lk max / a) 3 times larger than the 
standard deviation. Since the horizon size at the end of 
preheating is much larger than the shortest wavelength 
in the resonance band, 

the required amplitude is extremely large, and hence 
the probability of finding it is completely negligible. 

Hence, non-Gaussianity cannot affect the estimate of 
the PBH formation rate. 



V. THE EFFECT OF METRIC 
PERTURBATIONS 

Finally we briefly discuss the role of gravitational inter- 
actions (metric perturbations) which might have effects 
on PBH formation. We have neglected them throughout 
this paper but there is a good reason why one expects 
that this is a good approximation. The time scale for 
gravitational collapse is at most free fall time. Unless 
density perturbations significantly exceed O(l), this time 
scale is identical to the Hubble time scale. On the other 
hand parametric resonance undergoes with the time scale 
determined by the effective mass of the inflation, which 
is in general shorter than the Hubble time scale. More- 
over, in the expanding universe gravitational instability 
does not grow exponentially, while parametric resonance 
drives exponential growth of perturbations. Hence, we 
expect gravitational interactions play a subdominant role 
at the stage of preheating, although later on a longer time 
scale gravitational collapse may proceed in cases where 
the effective pressure happens to be very small (such as 
in the massive model). 

In the treatment neglecting the gravitational interac- 
tion, there arises another subtlety related to the gauge. 
We discussed the amplitude of density perturbations at 
the horizon size, but it is a gauge dependent quantity. 
Hence, strictly speaking, it is incorrect to quote the cri- 
terion for the PBH formation stated in terms of den- 
sity perturbations in comoving gauge. Moreover, it is 
more suitable to use the amplitude of metric perturba- 
tions rather than density contrast [iol Elj . 

In the present context, this is mainly because the 
power spectrum of density perturbation has a strong 
/c-dependence, oc fc 3 , which means that probability of 
PBH formation is very sensitive to the choice of the 
horizon size |26| . In contrast, since the metric pertur- 
bations well inside the horizon is characterized by the 
Newton potential, the power spectrum of metric per- 
turbations will be proportional to A; -1 . On the other 
hand, it was shown in |44j that the curvature pertur- 
bation on the constant energy density hypersurfaces £ 
(Bardeen parameter) behaves like oc k 3 on super horizon 
scales after prehea ting by using the separate universe ap- 
proach m ei m El m ■ 

Therefore, say, the curvature perturbation £ will have 
a peak near at the horizon size, and we will be able to 
obtain an unambiguous upper limit on the abundance of 
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PBHs produced by preheating. From the above consid- 
eration, including metric perturbation in the evolution 
of scalar fields is an interesting issue which we leave for 
future work. 



VI. SUMMARY AND DISCUSSION 

We have studied the formation of black holes during 
preheating after inflation. Preheating provides a chal- 
lenging framework to illucidate various complex physical 
processes such as non-equilibrium, non-perturbative field 
theory in expanding backgrounds. In addition, since pre- 
heating is generic in some regions of parameter space for 
many inflationary models, the issue of whether primor- 
dial black holes (PBH) are overproduced is an important 
one. 

To address this we have performed two and three di- 
mensional lattice simulations of several different inflaton 
potentials (conformal and massive) and parameter re- 
gions. These simulations automatically incorporate all 
non-linear effects such as backreaction and rescattering 
of fields. We found no evidence for over-production of 
PBH, in contrast to earlier expectations. In addition we 
found that although highly non-Gaussian on very small 
scale, the spectrum of density perturbations is effectively 
Gaussain on the horizon (i.e. Hubble) scale. 

Our results can be understood simply. For all cases, 
we found that when the amplitude of density perturba- 
tions at about the shortest wavelength in the resonant 
band becomes of order unity, the growth due to para- 
metric resonance terminates due to backreaction. At the 
end of preheating the final density spectrum is universal, 
with a blue power spectrum, oc k D (D = 2,3) on hori- 
zon scales, depending on whether the simulation was two 
(D = 2) or three (D = 3) dimensional. Since the peak 
of the spectrum is scales significantly shorter than the 
horizon scale, the amplitude of the density perturbation 
at the horizon scale extrapolated from the peak is typ- 
ically about an order of magnitude below the threshold 
for PBH overproduction. 

We gave an explanation of this universality in the slope 
of the final power spectrum on large scales as a result of 
loss of coherence due to parametric resonance. 

These results argue for the view that gencrically PBHs 
will not be produced so much as to violate astrophys- 
ical constraints, even in the case of strong preheating. 
This result also applies to the cases with other models 
and parameters regions if our interpretation of the power 
spectrum on large scales is universally correct. 

The result obtained in thispaper is different from the 
claim by Green and Malik |3J], where they found that 
PBHs are likely to be overproduced. They estimated 
the time when backreaction becomes significant as well 
as the time when the amplitude of density perturbations 
exceeds the threshold separately based on linear approx- 
imation. Comparison of these two times was used to give 
a criteria for PBH formation. However, for example, in 



Ref. [13| the time when the backreaction becomes signif- 
icant in the case with m = 10 _6 m p j and g = 10~ 3 is 
estimated to be ~ 90, which is slightly later than our 
numerical result (See Fig. |SJ). Since the growth of per- 
turbation amplitude is exponential, a small error in the 
back reaction time can lead to wrong conclusions. In cal- 
culating the abundance of produced PBHs, only 10 ~ 20 
percent error of backreaction time can lead to the oppo- 
site conclusion. By contrast, our conclusion is based on 
self-consistent simulations and rather robust qualitative 
observations. There is no delicate comparison of different 
time scales. 

In this paper, we considered standard slow roll infla- 
tion. In such cases, preheating occurs at rather high 
energies. Therefore the mass of PBHs formed in the 
present context is too small to avoid evaporation before 
the Big Bang Nucleosynthesis. Hence, those PBHs are 
not subject to any observational constraint even if PBHs 
are produced abundantly, unless PBHs leave Planck mass 
relics @. (Hence we have assumed that PBH will leave 
Planck mass relic throughout this paper. ) Even if relics 
are formed, the constraint on the mass fraction of pro- 
duced black holes is (3 < 10~ 20 0, we can therefore con- 
clude that production of PBHs by preheating does not 
give a serious constraint on such simple models of pre- 
heating as discussed in this paper. 

However, our qualitative results will also apply for pre- 
heating at lower energies. Let us consider the possibility 
of making more massive black holes by preheating, say, 
in the case of massive inflaton models. Now we consider 
to vary the inflaton mass m and the value of </>-filed at 
the beginning of preheating within 4>o ~ n^pi (As an 
example of realizing small value of cf>o, we can consider 
the hybrid inflation model 49].). The key quantity is 
the ratio of k max given in (|13fl to the Hubble parameter, 
which is estimated as 

kmax ITlpl 1/4 ,„„s 



This ratio is independent of m and is the smallest for 
<t> ~ ni.pi. Then the same estimate given in (13411 applies 
as an upper bound for 5h- On the other hand, lowering 
(3 down to 10~ 30 only change the upper bound on a from 
0.031 to 0.026. Thus we can say that overproduction of 
more massive PBHs due to parametric resonance is also 
unlikely. 

We have also confirmed Gaussianity of the probability 
distribution of density perturbations at the horizon scale, 
which is assumed in the estimate of the production rate 
of PBHs. The appearance of Gaussianity is in accordance 
with the interpretation of the spectrum on large scales. 
If perturbations become uncorrelated beyond the short- 
est resonant scale, perturbations at the horizon scale are 
given by average of many statistically independent ran- 
dom variables. Thus Gaussianity naturally follows from 
the central limit theorem. 
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